// Growth rate calculator.
// Code that calculates eigenvalue for LexminDFA.

#include "grate/eigenvalue.h"
#include "grate/scc_builder.h"

#include <cmath>
#include "grate/dfa.h"

using growth_rate::LexminDFA;
typedef growth_rate::DFAInterface::Edge DFAEdge;

namespace growth_rate {

static const size_t kMinIterations = 25;

static const double kDefaultPrecision = 1e-6;
static const size_t kDefaultMaxIterations = 20000;

EigenvalueOptions::EigenvalueOptions()
    : scc_strategy(FindSccDefault),
      precision(kDefaultPrecision),
      max_iterations(kDefaultMaxIterations),
      add_loops_to_nodes(false),
      use_oracle(false),
      verbosity(false),
      generate_edge_list_mode(EdgeListModeAuto) {}

// Generates edge list of a given DFA it is specified by options.
static bool get_edge_list(
    const DFAInterface *dfa, const EigenvalueOptions &options,
    ScopedArray<DFAEdge> *edge_list, size_t num_edges) {
  // Automatic detection tells that we need to generate the edge list if
  // it helps us to save memory. For this both conditions should hold:
  //   (1) edge list size is less than arrays of "double" counters used
  //       to calculate eigenvalues.
  //   (2) edge list size is less than transition table size.
  size_t edge_list_size = sizeof(DFAEdge) * num_edges;
  size_t counters_size = 2 * sizeof(double) * dfa->size();  // NOLINT
  if (options.generate_edge_list_mode == EigenvalueOptions::EdgeListModeAlways
      || (options.generate_edge_list_mode == EigenvalueOptions::EdgeListModeAuto
          && edge_list_size < dfa->transition_table_size()
          && edge_list_size < counters_size)) {
    edge_list->initialize(num_edges);
    dfa->dump_valid_edges(edge_list->data());
    return true;
  }
  return false;
}

double EigenvalueCalculator::get_max_eigenvalue(
    DFAInterface *dfa, const EigenvalueOptions &options,
    size_t *iterations_used) {
  // Delete DFA at exit from this function.
  AutoObjectDeleter<DFAInterface> dfa_deleter(dfa);
  if (iterations_used != 0)
    *iterations_used = 0;

  CHECK_LT(options.scc_strategy, EigenvalueOptions::FindSccCnt);
  if (options.scc_strategy == EigenvalueOptions::FindSccNever) {
    CHECK(options.use_oracle &&
          "FindSccNever strategy only works with oracle!");
  } else if (options.scc_strategy == EigenvalueOptions::FindSccDefault) {
    // Build strongly connected components using naive Tarjan algorithm.
    SccPartition *scc_partition = SccBuilder::build_scc_naive(dfa);
    AutoObjectDeleter<SccPartition> scc_deleter(scc_partition);
    int scc_number = scc_partition->scc_number();
    if (scc_number == 0)
      return 0.0;
    if (scc_number > 1) {
      Report("Unexpected: DFA has %d SCC\n", scc_number);
      UNIMPLEMENTED();
      return 0.0;
    }
    dfa->leave_single_scc(scc_partition);
  } else {
    UNIMPLEMENTED();
    return 0.0;
  }

  size_t num_nodes = dfa->size();
  size_t num_edges = dfa->number_of_valid_edges();
  if (options.scc_strategy != EigenvalueOptions::FindSccNever) {
    CHECK_LE(num_nodes, num_edges);
    if (options.verbosity) {
      Report("Number of nodes in SCC: %ld\n", num_nodes);
      Report("Number of edges in SCC: %ld\n", num_edges);
    }
    if (num_nodes == num_edges) {
      // The only SCC is a simple cycle.
      return 1.0;
    }
  }

  ScopedArray<DFAEdge> edge_list;
  bool have_edge_list = get_edge_list(dfa, options, &edge_list, num_edges);
  if (have_edge_list) {
    // Later in the algorithm, use edge list instead of transition table.
    // Clear the transition table now.
    dfa->shrink_size(0);
  }

  ScopedArray<double> cnt_old(num_nodes);
  ScopedArray<double> cnt_new(num_nodes);
  if (options.verbosity)
    Report("Eigenvalue counters allocated!\n");

  CHECK_GT(options.max_iterations, kMinIterations);
  ScopedArray<double> rates;
  if (options.use_oracle)
    rates.initialize(options.max_iterations);
  double prev_rate = 1.0;
  double last_local_min = 0.0;
  double last_local_max = 0.0;

  bool found = false;
  double result = 0.0;

  size_t iter;
  for (iter = 0; iter < options.max_iterations; iter++) {
    // Initialize counters.
    if (iter == 0) {
      if (options.use_oracle) {
        cnt_old.clear();
        cnt_old[LexminDFA::StartNode] = 1.0;
      } else {
        for (size_t node = 0; node < num_nodes; node++)
          cnt_old[node] = 1.0;
      }
    } else {
      memcpy(cnt_old.data(), cnt_new.data(), sizeof(cnt_old[0]) * num_nodes);
    }
    cnt_new.clear();

    // Update counters.
    if (have_edge_list) {
      for (size_t i = 0; i < num_edges; i++)
        cnt_new[edge_list[i].to] += cnt_old[edge_list[i].from];
    } else {
      for (size_t node = 0; node < dfa->size(); node++) {
        for (uint8_t letter = 0; letter < dfa->alphabet_size(); letter++) {
          int next = dfa->get_next(node, letter);
          if (next >= LexminDFA::StartNode)
            cnt_new[next] += cnt_old[node];
        }
      }
    }
    if (options.add_loops_to_nodes) {
      for (size_t node = 0; node < num_nodes; node++)
          cnt_new[node] += cnt_old[node];
    }

    if (options.use_oracle) {
      double total_rate = 0;
      for (size_t node = 0; node < num_nodes; node++) {
        total_rate += cnt_new[node];
      }
      if (total_rate == 0) {
        found = true;
        result = 0.0;
        if (options.verbosity)
          Report("Oracle: zero growth rate!\n");
        break;
      }
      rates[iter] = total_rate / prev_rate;
      prev_rate = total_rate;
    }

    if (iter < kMinIterations)
      continue;

    if (options.use_oracle) {
      if (rates[iter - 1] >= rates[iter - 2] &&
          rates[iter - 1] >= rates[iter])
        last_local_max = rates[iter - 1];
      if (rates[iter - 1] <= rates[iter - 2] &&
          rates[iter - 1] <= rates[iter])
        last_local_min = rates[iter - 1];

      if (fabs(rates[iter] - rates[iter - 1]) < options.precision) {
        if (last_local_max > 0 && last_local_min > 0) {
          // Oscillating sequence of rates.
          if (last_local_max - last_local_min < options.precision) {
            found = true;
            result = rates[iter];
            if (options.verbosity)
              Report("Oracle: oscillating sequence!\n");
            break;
          }
        } else {
          // Monotonic sequence of rates.
          found = true;
          double q = pow(fabs(rates[iter] - rates[iter - 1]), 1.0 / iter);
          double sum = (rates[iter] - rates[iter - 1]) / (1 - q);
          result = rates[iter - 1] + sum;
          if (options.verbosity)
            Report("Oracle: monotonic sequence!\n");
          break;
        }
      }
    } else {
      double minrate = 1e100;
      double maxrate = 0;
      for (size_t node = 0; node < num_nodes; node++) {
        double crate = cnt_new[node] / cnt_old[node];
        minrate = Min(minrate, crate);
        maxrate = Max(maxrate, crate);
      }

      if (maxrate - minrate < 2 * options.precision) {
        result = (minrate + maxrate) / 2.0;
        found = true;
        break;
      }
    }
  }

  if (found && options.add_loops_to_nodes)
    result--;

  if (iterations_used != 0)
    *iterations_used = iter + 1;
  return (found) ? result: -1.0;
}

}  // namespace growth_rate
